Proyecto: Riesgo en el Banco Giturra¶

MDS7202: Laboratorio de Programación Científica para Ciencia de Datos

Cuerpo Docente:¶

  • Profesor: Pablo Badilla, Ignacio Meza De La Jara
  • Auxiliar: Sebastián Tinoco
  • Ayudante: Diego Cortez M., Felipe Arias T.

Por favor, lean detalladamente las instrucciones de la tarea antes de empezar a escribir.

Equipo: SUPER IMPORTANTE¶

  • Nombre de alumno 1: Johnny Godoy
  • Nombre de alumno 2: Los terrores nocturnos de Johnny Godoy
  • Link de repositorio de GitHub: https://github.com/johnny-godoy/laboratorios-mds/blob/main/2023/proyecto2

Reglas¶

  • Fecha de entrega: 01/06/2021
  • Grupos de 2 personas.
  • Cualquier duda fuera del horario de clases al foro. Mensajes al equipo docente serán respondidos por este medio.
  • Estrictamente prohibida la copia.
  • Pueden usar cualquier material del curso que estimen conveniente.

Presentación del Problema¶

Giturra, un banquero astuto y ambicioso, estableció su propio banco con el objetivo de obtener enormes ganancias. Sin embargo, su reputación se vio empañada debido a las tasas de interés usureras que imponía a sus clientes. A medida que su banco crecía, Giturra enfrentaba una creciente cantidad de préstamos impagados, lo que amenazaba su negocio y su prestigio.

Para abordar este desafío, Giturra reconoció la necesidad de reducir los riesgos de préstamo y mejorar la calidad de los préstamos otorgados. Decidió aprovechar la ciencia de datos y el análisis de riesgo crediticio. Contrató a un equipo de expertos para desarrollar un modelo predictivo de riesgo crediticio.

Cabe señalar que lo modelos solicitados por el banquero deben ser interpretables. Ya que estos le permitira al equipo comprender y explicar cómo se toman las decisiones crediticias. Utilizando visualizaciones claras y explicaciones detalladas, pudieron identificar las características más relevantes, le permitirá analizar la distribución de la importancia de las variables y evaluar si los modelos son coherentes con el negocio.

Para esto Giturra les solicita crear un modelo de riesgo disponibilizandoles una amplia gama de variables de sus usuarios: como historiales de crédito, ingresos y otros factores financieros relevantes, para evaluar la probabilidad de incumplimiento de pago de los clientes. Con esta información, Giturra podra tomar decisiones más informadas en cuanto a los préstamos, ofreciendo condiciones más favorables a aquellos con menor riesgo de impago.

Instalación de Librerías y Carga de Datos.¶

Para el desarrollo de su proyecto, utilice el conjunto de datos dataset.pq para entrenar un modelo de su elección. Además, se adjunta junto con los datos del proyecto un archivo llamado requirements.txt que contiene todas las bibliotecas y versiones necesarias para el desarrollo del proyecto. Se le recomienda levantar un ambiente de conda para instalar estas librerías y así evitar cualquier problema con las versiones.


1. Introducción [0.5 puntos]¶

En este problema, se quiere clasificar si es un cliente mal candidato para un préstamo. En el conjunto de datos original, el problema es de 3 clases, pero debido al problema de negocios, se quiere resolver como un problema de clasificación binaria.

También se eliminó la información temporal de los datos, que originalmente se dividía por meses, y nos interesaba realizar predicciones a futuro sobre los clientes con información pasada, parece que esto se simplificó, reduciendo problemáticas de contaminación de datos y de sesgo temporal.

En particular, se quiere predecir si su puntaje es Poor o si no lo es, por lo que se plantea como un problema de clasificación binaria, a pesar de que tengamos más de dos categorías, aquellas que no son Poor se consideran como una nueva categoría Not Poor para este problema.

Las entradas de los datos incluyen información sobre el cliente (como su edad, sexo, y trabajo) sobre el préstamo (tasa de interés, duración, monto, etc.).

Debido a que se resuelve un problema de clasificación binaria desbalanceado, se utilizará el f1-score como métrica de evaluación de los modelos. Esto es porque no se tiene una métrica del "costo" de un error de clasificación que nos permita pesar los falsos positivos y falsos negativos de manera distinta. Sabemos que f1-score es robusto al desbalanceo de clases, y que es una métrica que combina precision y recall.

Finalmente se utilizó un modelo ExplainableBoostingMachine (EBM), que combina técnicas de ensamblaje de árboles de decisión (bagging y boosting) para obtener modelos potentes, incorporando restricciones que además hacen que el modelo sea mucho más interpretable. Este modelo fue el más preciso a pesar de ser tan interpretable como un modelo lineal. Las únicas transformaciones previas a los datos fueron una imputación simple por mediana en los datos numéricos faltantes y la creación de una nueva variable como el producto de dos otras, que fue altamente predictiva.

Este modelo logró resolver el problema existosamente, con 0.8 de f1 score final. Creo que lo único particularmente innovador que se realizó fue el uso de la variable producto. Otros equipos que hayan encontrado esta variable, u otras mejores, probablemente lograron mejores puntajes, dado que no se utilizaron técnicas de rebalanceo de datos, y se usó una grilla pequeña de datos. Otros modelos de boosting en general logran mejores desemepeños que EBM, así que si usaron una grilla más grande para un modelo como ImbalancedRandomForest, esperaría mejores resultados al costo de tener que usar metodologías de interpretación post-hoc, que son más limitadas a los gráficos que provee EBM.

2. Carga de datos Análisis Exploratorio de Datos [Sin puntaje]¶

Voy a escribir cosas en esta sección para poder recordar mejor.

In [1]:
import pickle
import warnings

import lightgbm as lgb
import numpy as np
import pandas as pd
import plotly.express as px
import seaborn as sns
import xgboost as xgb
from interpret import show
from interpret.glassbox import ExplainableBoostingClassifier
from optuna.distributions import IntDistribution
from optuna.integration.sklearn import OptunaSearchCV
from sklearn.compose import make_column_transformer
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import mutual_info_classif
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, roc_auc_score
from sklearn.model_selection import (
    GridSearchCV,
    StratifiedKFold,
    cross_validate,
    train_test_split,
)
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import Pipeline, make_pipeline
from sklearn.preprocessing import (
    FunctionTransformer,
    KBinsDiscretizer,
    OneHotEncoder,
    RobustScaler,
)
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.utils import resample
Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
In [2]:
df_full = pd.read_parquet("input/dataset.pq")
df_full.drop(columns=["customer_id"], inplace=True)
objects = df_full.select_dtypes("object")
object_columns = objects.columns
categories = df_full[object_columns].astype("category")
to_change = object_columns[
    categories.memory_usage(index=False) < objects.memory_usage(index=False)
]
df_full[to_change] = categories[to_change]
X, y = df_full.drop(columns=["credit_score"]), df_full["credit_score"]
In [3]:
df_full.describe()
Out[3]:
age annual_income monthly_inhand_salary num_bank_accounts num_credit_card interest_rate num_of_loan delay_from_due_date num_of_delayed_payment changed_credit_limit num_credit_inquiries outstanding_debt credit_utilization_ratio credit_history_age total_emi_per_month amount_invested_monthly monthly_balance credit_score
count 12500.000000 1.250000e+04 10584.000000 12500.000000 12500.000000 12500.000000 12500.000000 12500.000000 11660.00000 12246.000000 12243.000000 12500.000000 12500.000000 11380.000000 12500.000000 11914.000000 1.214500e+04 12500.000000
mean 105.771840 1.616206e+05 4186.634963 16.939920 23.172720 73.213360 3.099440 21.060880 32.93542 10.398582 26.292330 1426.220376 32.349265 18.230404 1488.394291 638.798715 -2.744614e+22 0.288160
std 664.502705 1.297842e+06 3173.690362 114.350815 132.005866 468.682227 65.105277 14.863091 237.43768 6.799253 181.821031 1155.169458 5.156815 8.302078 8561.449910 2049.195193 3.024684e+24 0.452924
min -500.000000 7.005930e+03 303.645417 -1.000000 0.000000 1.000000 -100.000000 -5.000000 -3.00000 -6.490000 0.000000 0.230000 20.100770 0.000000 0.000000 0.000000 -3.333333e+26 0.000000
25% 25.000000 1.945333e+04 1622.408646 3.000000 4.000000 8.000000 1.000000 10.000000 9.00000 5.370000 4.000000 566.072500 28.066517 12.000000 31.496968 73.736810 2.701501e+02 0.000000
50% 33.000000 3.757238e+04 3087.595000 6.000000 5.000000 14.000000 3.000000 18.000000 14.00000 9.410000 6.000000 1166.155000 32.418953 18.000000 72.887628 134.093193 3.393885e+02 0.000000
75% 42.000000 7.269021e+04 5967.937500 7.000000 7.000000 20.000000 5.000000 28.000000 18.00000 14.940000 10.000000 1945.962500 36.623650 25.000000 169.634826 261.664256 4.714245e+02 1.000000
max 8678.000000 2.383470e+07 15204.633333 1756.000000 1499.000000 5789.000000 1495.000000 67.000000 4293.00000 36.970000 2554.000000 4998.070000 48.199824 33.000000 81971.000000 10000.000000 1.463792e+03 1.000000
In [4]:
df_full.nunique()
Out[4]:
age                           258
occupation                     16
annual_income               12489
monthly_inhand_salary       10579
num_bank_accounts             174
num_credit_card               284
interest_rate                 300
num_of_loan                    77
delay_from_due_date            73
num_of_delayed_payment        129
changed_credit_limit         2985
num_credit_inquiries          201
outstanding_debt            12203
credit_utilization_ratio    12500
credit_history_age             34
payment_of_min_amount           3
total_emi_per_month         11226
amount_invested_monthly     11353
payment_behaviour               7
monthly_balance             12145
credit_score                    2
dtype: int64
In [5]:
df_full.dtypes
Out[5]:
age                          float64
occupation                  category
annual_income                float64
monthly_inhand_salary        float64
num_bank_accounts              int64
num_credit_card                int64
interest_rate                  int64
num_of_loan                  float64
delay_from_due_date            int64
num_of_delayed_payment       float64
changed_credit_limit         float64
num_credit_inquiries         float64
outstanding_debt             float64
credit_utilization_ratio     float64
credit_history_age           float64
payment_of_min_amount       category
total_emi_per_month          float64
amount_invested_monthly      float64
payment_behaviour           category
monthly_balance              float64
credit_score                   int64
dtype: object

df_numeric = df_full.select_dtypes(include=["float64", "int64"]) fig, axes = plt.subplots( nrows=len(df_numeric.columns), figsize=(10, 10*len(df_numeric.columns)) ) for col, ax in zip(df_numeric.columns, axes): try: sns.histplot( data=df_numeric, x=col, hue="credit_score", multiple="stack", ax=ax ) except ValueError: warnings.warn(f"Could not plot {col}")

  • Algunos histogramas muestran outliers enormes, que los vuelven ilegibles.

Claramente, algunos de estos casos son por imputación

  • No existe un histograma que trivialmente resuelva el problema.

Age Represents the age of the person

Occupation Represents the occupation of the person

Annual_Income Represents the annual income of the person

Monthly_Inhand_Salary Represents the monthly base salary of a person

Num_Bank_Accounts Represents the number of bank accounts a person holds

Interest_Rate Represents the interest rate on credit card

Num_of_Loan Represents the number of loans taken from the bank

Delay_from_due_date Represents the average number of days delayed from the payment date

Num_of_Delayed_Payment Represents the average number of payments delayed by a person

Changed_Credit_Limit Represents the percentage change in credit card limit

Num_Credit_Inquiries Represents the number of credit card inquiries

Credit_Mix Represents the classification of the mix of credits

Outstanding_Debt Represents the remaining debt to be paid (in USD)

Credit_History_Age Represents the age of credit history of the person

Payment_of_Min_Amount Represents whether only the minimum amount was paid by the person

Total_EMI_per_month Represents the monthly EMI payments (in USD)

Amount_invested_monthly Represents the monthly amount invested by the customer (in USD)

Payment_Behaviour Represents the payment behavior of the customer (in USD)

Monthly_Balance Represents the monthly balance amount of the customer (in USD)

3. Preparación de Datos [0.5 puntos]¶

3.1 Preprocesamiento con ColumnTransformer¶

Primero se crea una transformación de limpieza, que realiza codificaciones específicas que no se lograrían con un OrdinalEnconder sin definir un orden para las filas. Esta parte también se termina combinando algo con la ###

  • Se crea una nueva columna payment_behaviour_low_spent que codifica que payment_behaviour tiene Low_spent como substring
  • Análogamente, se crea payment_behaviour_value. Esta es ordinal: -1 codifica si está Small_value, 1 codifica si está Large_value, y se tiene un 0 en cualquier otro caso (pues no se puede saber)
  • También se crea una columna dummy payment_behaviour_na que codifica si payment_behaviour toma el valor desconocido (no interpretable) de !@9#%8
  • Se codifica payment_of_min_amount de manera ordinal. Si bien Yes y No tienen su codificación binaria obvia, hay que lidiar con el valor nulo NM, el que se codificó como -1 para ser distinto.
In [6]:
def cleanup_transform(X):
    X_ = X.copy()
    X_["payment_behaviour_low_spent"] = X_.payment_behaviour.apply(
        lambda x: "Low_spent" in x
    ).astype(int)
    X_["payment_behaviour_value"] = X_.payment_behaviour.apply(
        lambda x: -1 if "Small_value" in x else 1 if "High_value" in x else 0
    )
    X_["payment_behaviour_na"] = (X_.payment_behaviour == "!@9#%8").astype(int)
    X_["payment_of_min_amount"].replace(
        {"No": 0.0, "Yes": 1.0, "NM": -1.0}, inplace=True
    )
    X_["payment_of_min_amount"] = X_["payment_of_min_amount"].astype(float)
    return X_


cleanup = FunctionTransformer(cleanup_transform)

Además de las transformaciones mencionadas, el resto de las variables discretas (occupation y payment_behaviour) no tienen un orden asociado, así que se codifican con OneHotEncoder. Para el escalamiento, ya se observó la gran cantidad de outliers, así que se usa RobustScaler. Con esto se prueba la salida del transformer.

In [7]:
transformer = make_column_transformer(
    (
        OneHotEncoder(
            drop="first",
            sparse_output=False,
            handle_unknown="ignore",
            dtype=int,
        ),
        ["occupation", "payment_behaviour"],
    ),
    (
        RobustScaler(),  # No se usa StandardScaler porque es sensible a outliers
        X.select_dtypes(include=["float64", "int64"]).columns,
    ),
    remainder="passthrough",
)
transformer = make_pipeline(cleanup, transformer)
transformer.set_output(transform="pandas")
With transform="pandas", `func` should return a DataFrame to follow the set_output API.
Out[7]:
Pipeline(steps=[('functiontransformer',
                 FunctionTransformer(func=<function cleanup_transform at 0x00000290A9AF9F30>)),
                ('columntransformer',
                 ColumnTransformer(remainder='passthrough',
                                   transformers=[('onehotencoder',
                                                  OneHotEncoder(drop='first',
                                                                dtype=<class 'int'>,
                                                                handle_unknown='ignore',
                                                                sparse_output=False),
                                                  ['occupation',
                                                   'payment_behaviour']),
                                                 ('robustscaler'...
                                                  Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance'],
      dtype='object'))]))])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Pipeline(steps=[('functiontransformer',
                 FunctionTransformer(func=<function cleanup_transform at 0x00000290A9AF9F30>)),
                ('columntransformer',
                 ColumnTransformer(remainder='passthrough',
                                   transformers=[('onehotencoder',
                                                  OneHotEncoder(drop='first',
                                                                dtype=<class 'int'>,
                                                                handle_unknown='ignore',
                                                                sparse_output=False),
                                                  ['occupation',
                                                   'payment_behaviour']),
                                                 ('robustscaler'...
                                                  Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance'],
      dtype='object'))]))])
FunctionTransformer(func=<function cleanup_transform at 0x00000290A9AF9F30>)
ColumnTransformer(remainder='passthrough',
                  transformers=[('onehotencoder',
                                 OneHotEncoder(drop='first',
                                               dtype=<class 'int'>,
                                               handle_unknown='ignore',
                                               sparse_output=False),
                                 ['occupation', 'payment_behaviour']),
                                ('robustscaler', RobustScaler(),
                                 Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance'],
      dtype='object'))])
['occupation', 'payment_behaviour']
OneHotEncoder(drop='first', dtype=<class 'int'>, handle_unknown='ignore',
              sparse_output=False)
Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance'],
      dtype='object')
RobustScaler()
passthrough
In [8]:
X_transformed = transformer.fit_transform(X)
X_transformed.head()
Out[8]:
onehotencoder__occupation_Architect onehotencoder__occupation_Developer onehotencoder__occupation_Doctor onehotencoder__occupation_Engineer onehotencoder__occupation_Entrepreneur onehotencoder__occupation_Journalist onehotencoder__occupation_Lawyer onehotencoder__occupation_Manager onehotencoder__occupation_Mechanic onehotencoder__occupation_Media_Manager ... robustscaler__outstanding_debt robustscaler__credit_utilization_ratio robustscaler__credit_history_age robustscaler__total_emi_per_month robustscaler__amount_invested_monthly robustscaler__monthly_balance remainder__payment_of_min_amount remainder__payment_behaviour_low_spent remainder__payment_behaviour_value remainder__payment_behaviour_na
0 0 0 0 0 0 0 0 0 0 0 ... -0.258118 -0.991589 NaN -0.168764 -0.581650 0.093085 0.0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 ... -0.406645 0.060172 0.692308 -0.391431 0.451297 0.082920 0.0 1 -1 0
2 0 0 0 1 0 0 0 0 0 0 ... 0.099178 0.696004 0.000000 1.260369 52.498488 2.762925 0.0 0 -1 0
3 0 0 0 0 1 0 0 0 0 0 ... -0.386766 -0.594409 -0.076923 -0.408810 -0.045102 0.197878 0.0 0 -1 0
4 0 1 0 0 0 0 0 0 0 0 ... -0.161096 -0.766148 1.000000 -0.527644 0.251361 0.122278 1.0 0 -1 0

5 rows × 42 columns

Todo cuadra.

3.2 Holdout¶

Generaré conjuntos de validación de todos modos para elegir cuales modelos usar en la búsqueda de hiperparámetros, pues no tiene sentido realizar la elección con respecto al conjunto de prueba. Esto permitirá que toda llamada de cross_validate, GridSearchCV y OptunaGridSearchCV sea consistente por utilizar los mismos conjuntos de validación (en principio, debería serlo igual debido a que tienen el mismo parámetro por defecto en validación cruzada, pero no está mal asegurarse en caso de que se rompa eso en versiones futuras).

In [9]:
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0, stratify=y
)
kf = StratifiedKFold(shuffle=True, random_state=0)

Se realizará adversarial cross validation, es decir, se intentará predecir si un dato pertenece al conjunto de entrenamiento o al conjunto de validación. El puntaje auc debería ser cercano a 0.5 si el modelo no puede distinguir entre ambos conjuntos, que es lo deseado.

Si se tiene un puntaje más alto, existe un problema de data leakage. Para encontrar la razón, se usa un modelo interpretable ExplainableBoostingMachine.

In [10]:
ebmclf = ExplainableBoostingClassifier(random_state=0)
X_adv = pd.concat([X_train, X_test])
y_adv = pd.Series(
    np.concatenate([np.zeros(len(X_train)), np.ones(len(X_test))]),
    index=X_adv.index,
)
ebmclf.fit(X_adv, y_adv)
roc_auc_score(y_adv, ebmclf.predict_proba(X_adv)[:, 1])
Missing values detected. Our visualizations do not currently display missing values. To retain the glassbox nature of the model you need to either set the missing values to an extreme value like -1000 that will be visible on the graphs, or manually examine the missing value score in ebm.term_scores_[term_index][0]
Out[10]:
0.6059684

El puntaje es cercano a 0.5, por lo que la preocupación de data leakage es menor, pero no se puede descartar. Se procede a explorar las variables más importantes.

In [11]:
show(ebmclf.explain_global())

Los efectos encontrados no son preocupantes. En general los puntajes (que se interpretan igual a un valor de SHAP, como logits) son muy cercanos a cero, y las pocas secciones de los gráficos que no son cercanas a cero tienen muy alta incertidumbre al correr el modelo en distintos submuestreos. Se procede a usar el holdout propuesto.

También se verifican los conjuntos de validación. No se intenciona realiza un análisis más profundo, así que se usa LightGBM el puntaje auc de manera veloz y sencilla, mostrando la variable más importante, pues si el puntaje auc es alto, entonces esta variable puede estar contaminada.

In [12]:
df_adv_report = pd.DataFrame(
    {
        "auc": np.zeros(5),
        "most_important_feature": np.zeros(5, dtype="object"),
    },
    index=[f"Fold {i}" for i in range(5)],
)
for i, (train_index, test_index) in enumerate(kf.split(X_adv, y_adv)):
    X_train_cv, X_test_cv = X_adv.iloc[train_index], X_adv.iloc[test_index]
    y_train_cv, y_test_cv = y_adv.iloc[train_index], y_adv.iloc[test_index]
    lgbclf = lgb.LGBMClassifier(random_state=0)
    lgbclf.fit(
        X_train_cv,
        y_train_cv,
    )
    pred = lgbclf.predict_proba(X_test_cv)[:, 1]
    df_adv_report.loc[f"Fold {i}", "auc"] = roc_auc_score(y_test_cv, pred)
    df_adv_report.loc[f"Fold {i}", "most_important_feature"] = X_train_cv.columns[
        lgbclf.feature_importances_.argmax()
    ]
df_adv_report
Out[12]:
auc most_important_feature
Fold 0 0.482397 changed_credit_limit
Fold 1 0.489064 changed_credit_limit
Fold 2 0.515999 monthly_inhand_salary
Fold 3 0.502982 monthly_inhand_salary
Fold 4 0.518786 monthly_inhand_salary

La validación adversarial no muestra problemas de data leakage, pues todos los puntajes son esencialmente azarosos.

3.3 Datos nulos.¶

Podemos ver que las únicas columnas con datos nulos son las numéricas. Por lo tanto, se usará una estrategia de imputación simple por la mediana, y agregando una columna indicatriz de los valores faltantes para no perder información.

Se usa la mediana en vez de la media porque ya observamos que hay fuertes outliers.

In [13]:
X_transformed.isna().sum().sort_values(ascending=False)
Out[13]:
robustscaler__monthly_inhand_salary                                  1916
robustscaler__credit_history_age                                     1120
robustscaler__num_of_delayed_payment                                  840
robustscaler__amount_invested_monthly                                 586
robustscaler__monthly_balance                                         355
robustscaler__num_credit_inquiries                                    257
robustscaler__changed_credit_limit                                    254
onehotencoder__occupation_Architect                                     0
robustscaler__num_bank_accounts                                         0
robustscaler__num_credit_card                                           0
robustscaler__interest_rate                                             0
robustscaler__num_of_loan                                               0
robustscaler__delay_from_due_date                                       0
robustscaler__outstanding_debt                                          0
onehotencoder__occupation_Developer                                     0
robustscaler__credit_utilization_ratio                                  0
robustscaler__total_emi_per_month                                       0
remainder__payment_of_min_amount                                        0
remainder__payment_behaviour_low_spent                                  0
remainder__payment_behaviour_value                                      0
robustscaler__annual_income                                             0
robustscaler__age                                                       0
onehotencoder__payment_behaviour_Low_spent_Small_value_payments         0
onehotencoder__payment_behaviour_Low_spent_Medium_value_payments        0
onehotencoder__occupation_Doctor                                        0
onehotencoder__occupation_Engineer                                      0
onehotencoder__occupation_Entrepreneur                                  0
onehotencoder__occupation_Journalist                                    0
onehotencoder__occupation_Lawyer                                        0
onehotencoder__occupation_Manager                                       0
onehotencoder__occupation_Mechanic                                      0
onehotencoder__occupation_Media_Manager                                 0
onehotencoder__occupation_Musician                                      0
onehotencoder__occupation_Scientist                                     0
onehotencoder__occupation_Teacher                                       0
onehotencoder__occupation_Writer                                        0
onehotencoder__occupation________                                       0
onehotencoder__payment_behaviour_High_spent_Large_value_payments        0
onehotencoder__payment_behaviour_High_spent_Medium_value_payments       0
onehotencoder__payment_behaviour_High_spent_Small_value_payments        0
onehotencoder__payment_behaviour_Low_spent_Large_value_payments         0
remainder__payment_behaviour_na                                         0
dtype: int64
In [14]:
null_columns = X_transformed.columns[X_transformed.isna().any()]
In [15]:
imputer = make_column_transformer(
    (
        SimpleImputer(strategy="median", add_indicator=True),
        null_columns,
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
)

imputer.set_output(transform="pandas")
X_transformed_imputed = imputer.fit_transform(X_transformed)
X_transformed_imputed.head()
Out[15]:
robustscaler__monthly_inhand_salary robustscaler__num_of_delayed_payment robustscaler__changed_credit_limit robustscaler__num_credit_inquiries robustscaler__credit_history_age robustscaler__amount_invested_monthly robustscaler__monthly_balance missingindicator_robustscaler__monthly_inhand_salary missingindicator_robustscaler__num_of_delayed_payment missingindicator_robustscaler__changed_credit_limit ... robustscaler__interest_rate robustscaler__num_of_loan robustscaler__delay_from_due_date robustscaler__outstanding_debt robustscaler__credit_utilization_ratio robustscaler__total_emi_per_month remainder__payment_of_min_amount remainder__payment_behaviour_low_spent remainder__payment_behaviour_value remainder__payment_behaviour_na
0 -0.290586 -0.888889 0.194357 -0.333333 0.000000 -0.581650 0.093085 0.0 0.0 0.0 ... -0.916667 0.25 -0.833333 -0.258118 -0.991589 -0.168764 0.0 0 0 0
1 -0.011416 -1.111111 -0.416928 -0.666667 0.692308 0.451297 0.082920 0.0 0.0 0.0 ... -0.666667 -0.50 -0.833333 -0.406645 0.060172 -0.391431 0.0 1 -1 0
2 2.094020 -0.888889 -0.241379 -0.500000 0.000000 52.498488 2.762925 0.0 0.0 0.0 ... -0.500000 0.00 -0.555556 0.099178 0.696004 1.260369 0.0 0 -1 0
3 -0.109332 -0.555556 -0.775340 -0.333333 -0.076923 -0.045102 0.197878 0.0 0.0 0.0 ... -0.833333 -25.75 -0.777778 -0.386766 -0.594409 -0.408810 0.0 0 -1 0
4 -0.053914 0.111111 -0.713689 -0.333333 1.000000 0.251361 0.122278 0.0 0.0 0.0 ... -0.750000 -25.75 -0.944444 -0.161096 -0.766148 -0.527644 1.0 0 -1 0

5 rows × 49 columns

3.4 Feature Engineering [Bonus - 0.5 puntos]¶

En esta sección es interesante seguir con el EDA para encontrar transformaciones con datos más limpios. Veremos la información mutua de cada variable con respecto a la variable objetivo para tener una idea de buenas transformaciones.

In [16]:
discrete = X_transformed_imputed.dtypes == int
mi = mutual_info_classif(X_transformed_imputed, y, discrete_features=discrete)
mi = pd.Series(mi, index=X_transformed_imputed.columns).sort_values()
mi.plot.barh(figsize=(10, 20));

Muchas de las variables, particularmente algunas one-hot, tienen información mutua nula, es decir, son univariablemente independientes de la variable objetivo. Esto inspira a concentrarse más en variables numéricas que en otras.

Podríamos graficar las variables de mayor información mutua. Sin embargo, estas podrían ser altamente correlacionadas. Para determinar esto, se verá la matriz de correlación de las 5 variables con mayor información mutua.

In [17]:
mi_top_five = mi.sort_values(ascending=False).head(5).index
mi_simple_names = mi_top_five.str.split("__").str[1]
X_simple_names = X_transformed_imputed.rename(columns=lambda x: x.split("__")[1])
corr = X_simple_names[mi_simple_names].corr()
px.imshow(corr)

No hay correlación mayor por la cual preocuparse. Por lo tanto, se procede a graficar la distribución conjunta de las dos variables con mayor MI coloreadas por el objetivo, a buscar ideas de variables a agregar.

In [18]:
sns.jointplot(
    data=X_transformed_imputed,
    x="robustscaler__outstanding_debt",
    y="robustscaler__interest_rate",
    hue=y_train,
);

No se vé un patrón claro relacionado por una fórmula simple, como la razón entre las variables. Se nota que cuando la deuda (escalada) es negativa, las tasas de interés son más concentradas, y hay menos clientes riegosos. Cuando es positiva, la mayoría de clientes con tasa de interés positiva son riesgosos. Por lo tanto, se agrega el producto de las variables, que tendrá el signo de la deuda (dado que la tasa queda casi siempre positiva) con una magnitud escalada por la tasa de interés.

In [19]:
X_train["debt_x_interest_rate"] = X_train["outstanding_debt"] * X_train["interest_rate"]
X_test["debt_x_interest_rate"] = X_test["outstanding_debt"] * X_test["interest_rate"]
X_transformed = transformer.fit_transform(X_train)
X_transformed_imputed = imputer.fit_transform(X_transformed)

Como se usarán ensamblajes en árboles, una buena transformación a realizar es discretización de las variables continuas. Esto regulariza a cada árbol y hace la búsqueda de splits mucho más veloz.

XGBoost y LightGBM la realizan por debajo, pero RandomForest no, por lo que se utilizará también KBinsDicretizer en ese caso.

In [20]:
discretizer = make_column_transformer(
    (
        KBinsDiscretizer(
            n_bins=256, encode="ordinal", random_state=0, dtype=np.float32
        ),
        X_transformed_imputed.columns[X_transformed_imputed.dtypes != int],
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
)
discretizer.set_output(transform="pandas")
discretizer.fit_transform(X_transformed_imputed).head()
C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 0 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 1 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 2 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 3 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 4 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 5 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 6 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 7 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 8 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 9 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 10 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 11 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 12 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 13 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 14 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 16 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 17 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 18 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 19 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 20 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 23 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 24 are removed. Consider decreasing the number of bins.

C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_discretization.py:313: UserWarning:

Bins whose width are too small (i.e., <= 1e-8) in feature 26 are removed. Consider decreasing the number of bins.

Out[20]:
robustscaler__monthly_inhand_salary robustscaler__num_of_delayed_payment robustscaler__changed_credit_limit robustscaler__num_credit_inquiries robustscaler__credit_history_age robustscaler__amount_invested_monthly robustscaler__monthly_balance missingindicator_robustscaler__monthly_inhand_salary missingindicator_robustscaler__num_of_delayed_payment missingindicator_robustscaler__changed_credit_limit ... onehotencoder__occupation_Writer onehotencoder__occupation________ onehotencoder__payment_behaviour_High_spent_Large_value_payments onehotencoder__payment_behaviour_High_spent_Medium_value_payments onehotencoder__payment_behaviour_High_spent_Small_value_payments onehotencoder__payment_behaviour_Low_spent_Large_value_payments onehotencoder__payment_behaviour_Low_spent_Medium_value_payments onehotencoder__payment_behaviour_Low_spent_Small_value_payments remainder__payment_behaviour_low_spent remainder__payment_behaviour_na
8746 3.0 27.0 240.0 6.0 1.0 0.0 36.0 0.0 0.0 0.0 ... 0 0 0 1 0 0 0 0 0 0
9550 44.0 17.0 210.0 12.0 9.0 115.0 51.0 0.0 1.0 0.0 ... 0 0 0 0 0 0 0 1 1 0
6961 214.0 5.0 47.0 6.0 31.0 228.0 234.0 0.0 0.0 0.0 ... 0 0 0 0 0 0 0 0 0 1
10085 45.0 30.0 122.0 12.0 15.0 82.0 70.0 0.0 0.0 0.0 ... 0 0 0 0 0 0 1 0 1 0
2864 141.0 17.0 24.0 10.0 18.0 66.0 190.0 0.0 1.0 0.0 ... 0 0 0 1 0 0 0 0 0 0

5 rows × 50 columns

Más áun, se intentará ver la distribución de las dos variables con mayor información mútua, coloreada por la variable objetivo.

4. Baseline [1.5 puntos]¶

In [21]:
transformer_with_imputer = make_pipeline(transformer, imputer)
transformer_with_imputer.set_output(transform="pandas")

dummy = DummyClassifier(strategy="stratified", random_state=0)
logreg = make_pipeline(
    transformer_with_imputer,
    LogisticRegression(
        random_state=0, class_weight="balanced", n_jobs=-2, solver="sag"
    ),
)
knn = make_pipeline(
    transformer_with_imputer, KNeighborsClassifier(n_jobs=-2, algorithm="kd_tree")
)
dt = make_pipeline(
    transformer,  # no necesita imputar
    DecisionTreeClassifier(random_state=0, class_weight="balanced"),
)
svc = make_pipeline(
    transformer_with_imputer,
    SVC(random_state=0, class_weight="balanced", probability=True),
)
rf = make_pipeline(
    transformer_with_imputer,
    discretizer,
    RandomForestClassifier(random_state=0, n_jobs=-2, class_weight="balanced"),
)
# No necesitan pipeline:
# realizan codificación categórica, lidian con valores nulos y no necesitan escalamiento
lgbclf = lgb.LGBMClassifier(random_state=0, n_jobs=-2)
xgbclf = xgb.XGBClassifier(
    random_state=0, enable_categorical=True, n_jobs=-2, tree_method="hist"
)
# Modelo extra a probar: EBM
simple_imputer = make_column_transformer(
    (SimpleImputer(strategy="median"), X_train.select_dtypes("number").columns),
    remainder="passthrough",
)
ebmclf = make_pipeline(
    simple_imputer,
    ExplainableBoostingClassifier(random_state=0),
)
models = {
    "Dummy": dummy,
    "LogisticRegression": logreg,
    "KNeighborsClassifier": knn,
    "DecisionTreeClassifier": dt,
    "SVC": svc,
    "RandomForestClassifier": rf,
    "LightGBMClassifier": lgbclf,
    "XGBClassifier": xgbclf,
    "ExplainableBoostingClassifier": ebmclf,
}
C:\Users\David\PycharmProjects\laboratorio-mds\venv\lib\site-packages\sklearn\preprocessing\_function_transformer.py:345: UserWarning:

With transform="pandas", `func` should return a DataFrame to follow the set_output API.

No encuentro que sea correcto utilizar el conjunto de test todavía, por lo que no usaré classification_report, sino que usaré cross_validate para obtener las métricas más robustas de validación. Ver el conjunto de test arruina el punto de usar GridSeachCV en la siguiente sección. Se reportarán las mismas métricas de classification_report para cada modelo.

In [22]:
report = {}
metrics = ["accuracy", "precision", "recall", "f1"]
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for name, model in models.items():
        print(name, f"\n{'-' * (len(name))}")
        cv = cross_validate(model, X_train, y_train, scoring=metrics, cv=kf)
        report[name] = cv
        print(pd.DataFrame(cv).mean(), "\n")
Dummy 
-----
fit_time          0.003600
score_time        0.039809
test_accuracy     0.595300
test_precision    0.304727
test_recall       0.315406
test_f1           0.309974
dtype: float64 

LogisticRegression 
------------------
fit_time          0.714425
score_time        0.040705
test_accuracy     0.503900
test_precision    0.339810
test_recall       0.710666
test_f1           0.452906
dtype: float64 

KNeighborsClassifier 
--------------------
fit_time          0.238784
score_time        2.122591
test_accuracy     0.736400
test_precision    0.555156
test_recall       0.430942
test_f1           0.485148
dtype: float64 

DecisionTreeClassifier 
----------------------
fit_time          0.838565
score_time        0.035408
test_accuracy     0.719400
test_precision    0.514518
test_recall       0.559000
test_f1           0.534165
dtype: float64 

SVC 
---
fit_time          28.748667
score_time         2.527197
test_accuracy      0.721400
test_precision     0.116783
test_recall        0.115972
test_f1            0.116376
dtype: float64 

RandomForestClassifier 
----------------------
fit_time          0.957615
score_time        0.090821
test_accuracy     0.782900
test_precision    0.671177
test_recall       0.485069
test_f1           0.562645
dtype: float64 

LightGBMClassifier 
------------------
fit_time          0.164837
score_time        0.018004
test_accuracy     0.788100
test_precision    0.666916
test_recall       0.530170
test_f1           0.590366
dtype: float64 

XGBClassifier 
-------------
fit_time          0.223851
score_time        0.016004
test_accuracy     0.775500
test_precision    0.637509
test_recall       0.512824
test_f1           0.568075
dtype: float64 

ExplainableBoostingClassifier 
-----------------------------
fit_time          4.207654
score_time        0.023054
test_accuracy     0.790500
test_precision    0.674245
test_recall       0.528782
test_f1           0.592309
dtype: float64 

Estudiemos específicamente la métrica F1, que es la que importará para la decisión final.

In [23]:
f1_vals = pd.DataFrame({name: cv["test_f1"] for name, cv in report.items()})
with open("output/f1_scores_baseline.pkl", "wb") as f:
    pickle.dump(f1_vals, f)
f1_vals.mean().sort_values(ascending=False)
Out[23]:
ExplainableBoostingClassifier    0.592309
LightGBMClassifier               0.590366
XGBClassifier                    0.568075
RandomForestClassifier           0.562645
DecisionTreeClassifier           0.534165
KNeighborsClassifier             0.485148
LogisticRegression               0.452906
Dummy                            0.309974
SVC                              0.116376
dtype: float64

También es interesante ver como se distribuyen los valores de F1 para cada modelo.

In [24]:
order_index = f1_vals.mean().sort_values().index
fig = f1_vals[order_index].boxplot(
    figsize=(20, 10), rot=45, showmeans=True, fontsize=20
)
fig.set_title("F1 scores for each model")
fig.set_ylabel("F1 score");
  • ¿Hay algún clasificador entrenado mejor que el azar (Dummy)? R: Todos los basados

en árboles de decisión, y el logístico. El único peor es SVC.

  • ¿Cuál es el mejor clasificador entrenado? R: ExplainableBoostingMachine, que tiene

alto desemepeño, y como se verá adelante, es fácil de interpretar.

  • ¿Por qué el mejor clasificador es mejor que los otros? R: Los clasificadores

basados en árboles pueden aprender relaciones no-lineales de manera no paramétrica sin sufrir la maldición de la dimensionalidad, a diferencia del resto de los modelos. Sin embargo, CART es un algoritmo ávaro que no encuentra un árbol óptimo en general. Variantes pequeñas tienen alto sesgo y variantes grandes tienen alta varianza. Los modelos de ensamblaje boosting y bagging resuelven esos problemas, respectivamente, por lo que logran mejores puntajes. Los dos mejores algoritmos utilizan boosting. El mejor algoritmo, ExplainableBoostingMachine tiene una forma funcional más restringida, evitando sobreajustes, y además utiliza bagging para reducir su varianza.

  • Respecto al tiempo de entrenamiento, con cual cree que sería mejor experimentar R:

LightGBM es un modelo de boosting con una implementación extremadamente veloz, así que buscar hiperparámetros no es muy problemático. En mi experiencia, he visto sin embargo que es menos sensible a cambios de hiperparámetros que XGBoost, pero sí más lento. ExplainableBoostingMachine desafortunadamente es bastante lento de entrenar, pero su guía de hiperparámetros es sencilla y dice que debe cambiar muy pocos parámetros.

5. Optimización del Modelo [1.5 puntos]¶

Para crear las grillas, se usan las guías de hiperparámetros de cada uno delos modelos. LGB usa una grilla más grande debido a su alta velocidad de entrenamiento, que está basada en su guía de hiperparámetros.

In [25]:
ebm_grid = {"max_leaves": [2, 3, 5], "max_bins": [32, 256, 1024]}
lgb_grid = {
    "extra_trees": [True, False],
    "max_bin": [127, 255, 511],
    "num_leaves": [31, 63, 127],
    "boosting_type": ["gbdt", "dart", "goss"],
}
In [26]:
def get_gridsearch(pipeline_or_model, params):
    if isinstance(pipeline_or_model, Pipeline):
        model_name = pipeline_or_model.steps[-1][0]
        model_params = {f"{model_name}__{k}": v for k, v in params.items()}
    else:
        model_params = params
    return GridSearchCV(
        pipeline_or_model,
        model_params,
        cv=kf,
        scoring="f1",
    )
In [27]:
ebm_best = get_gridsearch(
    ebmclf,
    ebm_grid,
)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    ebm_best.fit(X_train, y_train)
In [28]:
lgb_best = get_gridsearch(
    lgbclf,
    lgb_grid,
)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    lgb_best.fit(X_train, y_train)
In [29]:
scores = pd.DataFrame(
    {
        "Model": ["EBM", "LGBM"],
        "F1": [ebm_best.best_score_, lgb_best.best_score_],
    }
)
scores
Out[29]:
Model F1
0 EBM 0.595611
1 LGBM 0.599526

La diferencia no fue enorme, con ambos modelos cambiando su puntaje solo marginalmente. LGBM se acercó más a 0.6, invirtiendo el orden de los modelos. Sin embargo, no es una diferencia significativa como para justificar la pérdida de interpretabilidad. Se realizará entonces una búsqueda con Optuna en EBM, a ver si esto cambia los resultados.

In [30]:
ebm_opt_grid = {
    "explainableboostingclassifier__max_leaves": IntDistribution(2, 5),
    "explainableboostingclassifier__max_bins": IntDistribution(32, 1024, log=True),
}
ebm_best_opt = OptunaSearchCV(
    ebmclf,
    ebm_opt_grid,
    cv=kf,
    scoring="f1",
    n_trials=15,
    timeout=600,
)
ebm_best_opt.fit(X_train, y_train)
C:\Users\David\AppData\Local\Temp\ipykernel_6388\1677874795.py:7: ExperimentalWarning:

OptunaSearchCV is experimental (supported from v0.17.0). The interface can change in the future.

[I 2023-07-15 23:01:01,026] A new study created in memory with name: no-name-c4a8609f-f8e9-4b94-872e-e16801b6e692
[I 2023-07-15 23:01:21,093] Trial 0 finished with value: 0.5916886599384277 and parameters: {'explainableboostingclassifier__max_leaves': 5, 'explainableboostingclassifier__max_bins': 463}. Best is trial 0 with value: 0.5916886599384277.
[I 2023-07-15 23:01:42,132] Trial 1 finished with value: 0.5886683739351001 and parameters: {'explainableboostingclassifier__max_leaves': 4, 'explainableboostingclassifier__max_bins': 829}. Best is trial 0 with value: 0.5916886599384277.
[I 2023-07-15 23:02:07,622] Trial 2 finished with value: 0.5944109344651529 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 111}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:02:27,526] Trial 3 finished with value: 0.5889441694608142 and parameters: {'explainableboostingclassifier__max_leaves': 4, 'explainableboostingclassifier__max_bins': 135}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:02:53,483] Trial 4 finished with value: 0.5922303867407761 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 189}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:03:18,230] Trial 5 finished with value: 0.5903504555454792 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 50}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:03:43,321] Trial 6 finished with value: 0.5931577156385999 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 195}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:04:09,342] Trial 7 finished with value: 0.5941240238338905 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 141}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:04:30,867] Trial 8 finished with value: 0.5929208231581513 and parameters: {'explainableboostingclassifier__max_leaves': 3, 'explainableboostingclassifier__max_bins': 375}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:04:57,732] Trial 9 finished with value: 0.5942900638602737 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 437}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:05:21,406] Trial 10 finished with value: 0.5870013637450819 and parameters: {'explainableboostingclassifier__max_leaves': 3, 'explainableboostingclassifier__max_bins': 41}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:05:44,413] Trial 11 finished with value: 0.5912835787261586 and parameters: {'explainableboostingclassifier__max_leaves': 3, 'explainableboostingclassifier__max_bins': 78}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:06:10,843] Trial 12 finished with value: 0.5932566953552102 and parameters: {'explainableboostingclassifier__max_leaves': 2, 'explainableboostingclassifier__max_bins': 315}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:06:32,412] Trial 13 finished with value: 0.5891662091065623 and parameters: {'explainableboostingclassifier__max_leaves': 3, 'explainableboostingclassifier__max_bins': 85}. Best is trial 2 with value: 0.5944109344651529.
[I 2023-07-15 23:06:54,560] Trial 14 finished with value: 0.590822376299512 and parameters: {'explainableboostingclassifier__max_leaves': 5, 'explainableboostingclassifier__max_bins': 1022}. Best is trial 2 with value: 0.5944109344651529.
Out[30]:
OptunaSearchCV(cv=StratifiedKFold(n_splits=5, random_state=0, shuffle=True),
               estimator=Pipeline(steps=[('columntransformer',
                                          ColumnTransformer(remainder='passthrough',
                                                            transformers=[('simpleimputer',
                                                                           SimpleImputer(strategy='median'),
                                                                           Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'd...
      dtype='object'))])),
                                         ('explainableboostingclassifier',
                                          ExplainableBoostingClassifier(random_state=0))]),
               n_jobs=1, n_trials=15,
               param_distributions={'explainableboostingclassifier__max_bins': IntDistribution(high=1024, log=True, low=32, step=1),
                                    'explainableboostingclassifier__max_leaves': IntDistribution(high=5, log=False, low=2, step=1)},
               scoring='f1', timeout=600)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
OptunaSearchCV(cv=StratifiedKFold(n_splits=5, random_state=0, shuffle=True),
               estimator=Pipeline(steps=[('columntransformer',
                                          ColumnTransformer(remainder='passthrough',
                                                            transformers=[('simpleimputer',
                                                                           SimpleImputer(strategy='median'),
                                                                           Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'd...
      dtype='object'))])),
                                         ('explainableboostingclassifier',
                                          ExplainableBoostingClassifier(random_state=0))]),
               n_jobs=1, n_trials=15,
               param_distributions={'explainableboostingclassifier__max_bins': IntDistribution(high=1024, log=True, low=32, step=1),
                                    'explainableboostingclassifier__max_leaves': IntDistribution(high=5, log=False, low=2, step=1)},
               scoring='f1', timeout=600)
Pipeline(steps=[('columntransformer',
                 ColumnTransformer(remainder='passthrough',
                                   transformers=[('simpleimputer',
                                                  SimpleImputer(strategy='median'),
                                                  Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance', 'debt_x_interest_rate'],
      dtype='object'))])),
                ('explainableboostingclassifier',
                 ExplainableBoostingClassifier(random_state=0))])
ColumnTransformer(remainder='passthrough',
                  transformers=[('simpleimputer',
                                 SimpleImputer(strategy='median'),
                                 Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance', 'debt_x_interest_rate'],
      dtype='object'))])
Index(['age', 'annual_income', 'monthly_inhand_salary', 'num_bank_accounts',
       'num_credit_card', 'interest_rate', 'num_of_loan',
       'delay_from_due_date', 'num_of_delayed_payment', 'changed_credit_limit',
       'num_credit_inquiries', 'outstanding_debt', 'credit_utilization_ratio',
       'credit_history_age', 'total_emi_per_month', 'amount_invested_monthly',
       'monthly_balance', 'debt_x_interest_rate'],
      dtype='object')
SimpleImputer(strategy='median')
passthrough
ExplainableBoostingClassifier(random_state=0)
In [31]:
with open("output/ebm_opt.pkl", "wb") as f:
    pickle.dump(ebm_best_opt, f)
In [32]:
ebm_best_opt.best_score_
Out[32]:
0.5944109344651529

El puntaje empeoró con respecto a la búsqueda anterior. Con los hiperparámetros fijos, se evalúan los modelos en el conjunto de prueba.

In [34]:
ebm_pred = ebm_best.predict(X_test)
print(classification_report(y_test, ebm_pred))
              precision    recall  f1-score   support

           0       0.83      0.90      0.86      1780
           1       0.69      0.54      0.60       720

    accuracy                           0.80      2500
   macro avg       0.76      0.72      0.73      2500
weighted avg       0.79      0.80      0.79      2500

In [35]:
lgb_pred = lgb_best.predict(X_test)
print(classification_report(y_test, lgb_pred))
              precision    recall  f1-score   support

           0       0.82      0.90      0.86      1780
           1       0.68      0.52      0.59       720

    accuracy                           0.79      2500
   macro avg       0.75      0.71      0.72      2500
weighted avg       0.78      0.79      0.78      2500

Acá vemos un resultado feliz: Si bien EBM tuvo peor puntaje en validación, fueron mejores en el conjunto de prueba en todas las métricas, logrando un accuracy del 0.8, f1-score de 0.73, precision y recall de 0.76 y 0.72. Entonces, este modelo tuvo los mejores puntajes y será el más interpretable.

6. Interpretabilidad [1.0 puntos]¶

Para esta sección, es importante saber que una ExplainableBoostingMachine es un modelo aditivo, es decir, su predicción es de la forma: $$ g(\mathbb{E}[y]) = \beta_0 + \sum_{i=1}^n f_i(x_i) + \sum_{i, j\in S} f_{ij}(x_i, x_j) $$

Donde $S$ es un conjunto de pares de características aprendido, y $g$ es la función logit en el caso de clasificación.

La capacidad de interpretación viene de que cada función $f_i$ y $f_{ij}$ se puede graficar, dando una descripción gráfica exacta del modelo sin técnicas de interpretación post-hoc:

  • El gráfico de $f_i$ da una interpretación de cómo

funciona el modelo marginalmente para la variable $i$. Esto es más informativo que un PDP, pues no es una marginalización, sino que da el efecto exacto que tiene un valor $x_i$ en el modelo.

  • El vector de los valores de SHAP para explicar un dato $x$ es precisamente $(f_i

(x_i))$, así que se puede calcular de manera exacta, sin usar aproximaciones.

  • La importancia total de la característica $i$ también se puede calcular fácilmente

como la media absoluta de $f_i$. Esto no depende del puntaje del modelo como la importancia de permutación, no cae en los sesgos de importancia basada en impureza, y es consistente con los valores de SHAP.

  • Como el modelo usa bagging, los puntajes $f(x_i)$ se calculan varias veces, por lo que se pueden calcular bandas de error que nos dan una idea de la incertidumbre de la explicación (que se traduce directamente en una incertidumbre del modelo).

Sin embargo, estas visualizaciones no funcionan muy bien para los Pipelines, de forma que se realizará de nuevo la imputación sin pipeline para esta sección.

In [36]:
simple_imputer = make_column_transformer(
    (SimpleImputer(strategy="median"), X_train.select_dtypes("number").columns),
    remainder="passthrough",
    verbose_feature_names_out=False,
)
simple_imputer.set_output(transform="pandas")
X_train_imputed = simple_imputer.fit_transform(X_train)
X_test_imputed = simple_imputer.transform(X_test)
ebm = ExplainableBoostingClassifier(random_state=0, max_leaves=2)
ebm.fit(X_train_imputed, y_train)
Out[36]:
ExplainableBoostingClassifier(max_leaves=2, random_state=0)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
ExplainableBoostingClassifier(max_leaves=2, random_state=0)
In [37]:
explanation = ebm.explain_global(name="EBM")
show(explanation)

Las variables más importantes fueron las que se vieron en el EDA con mayor información mutua, y la variable que se creó en la ingeniería de características multiplicando la deuda con la tasa de interés. En general, todas las características importantes son lógicas, que muestran datos de créditos y pagos del cliente, lo que se esperaría. El gráfico por defecto solamente muestra las más importantes, y ahí ya se nota que la distribución de importancia no es equitativa, lo que tiene sentido. El resto de variables es aún menos importante.

In [38]:
sample = resample(
    X_test_imputed, n_samples=10, random_state=0, stratify=y_test, replace=False
)
show(ebm.explain_local(sample, name="EBM"))

Los gráficos para cada observación muestran que el modelo es coherente, pues las variables más importantes son también siempre las monetarias. En el primer caso, primero se nota que el "valor base" del modelo es muy negativo, pues hay muy pocos casos de clientes riesgosos. Sin embargo, este cliente tiene alta deuda y tasa de interés (y aún más alta variable de producto que se creó) de 2300 y 53800 respectivamente, aumentando su riesgo altamente. Muy pocos factores siquiera reducen esta chace, como su changed_credit_limit.

En contraste, el segundo cliente tiene una mucha menor chance de ser riesgoso, con una deuda de 439 e interés de 1, factores que reducen su riesgo. Tiene factores que aumentan su riesgo, como el hecho de tener 2 pagos atrasados (que puede ser mal visto), pero no es suficientemente alto como para que sea riesgoso dado el valor base.

Para ver si hay variables irrelevantes, se imprimen todas las importancias.

In [39]:
importances = pd.Series(ebm.term_importances(), index=ebm.term_names_)
importances.sort_values()
Out[39]:
monthly_inhand_salary                           0.018784
occupation                                      0.018906
annual_income                                   0.022391
credit_utilization_ratio                        0.023124
age                                             0.024516
changed_credit_limit & debt_x_interest_rate     0.026860
delay_from_due_date & outstanding_debt          0.027019
interest_rate & debt_x_interest_rate            0.027712
interest_rate & delay_from_due_date             0.028924
delay_from_due_date & changed_credit_limit      0.034103
num_bank_accounts & delay_from_due_date         0.040404
delay_from_due_date & num_of_delayed_payment    0.040433
payment_of_min_amount                           0.041907
changed_credit_limit & outstanding_debt         0.049972
num_bank_accounts                               0.054201
num_credit_card & outstanding_debt              0.056775
total_emi_per_month                             0.066131
credit_history_age & debt_x_interest_rate       0.072739
num_of_loan                                     0.078135
credit_history_age                              0.096803
num_credit_inquiries                            0.096958
monthly_balance                                 0.100939
amount_invested_monthly                         0.113364
num_of_delayed_payment                          0.122184
interest_rate                                   0.153071
num_credit_card                                 0.166341
payment_behaviour                               0.182987
changed_credit_limit                            0.183356
delay_from_due_date                             0.215908
outstanding_debt                                0.277563
debt_x_interest_rate                            0.335878
dtype: float64

Variables como la edad, la ocupación, salario y ganancia anual tienen muy baja importancia, pero no es nula. Esto es debido a que EBM no realiza selección de características (fun fact: mi tesis de sobre desarrollar una variante que sí lo hace, pero no estaba lista para este proyecto :D). Los efectos de interacciones también son bajos, pero no nulos, es decir, es posible que haya valor en simplificar más aún el modelo a uno aditivo de la forma:

$$ g(\mathbb{E}[y]) = \beta_0 + \sum_{i=1}^n f_i(x_i) $$

Que se puede hacer con el parámetro interactions=0.

Esto es algo bueno, pues el modelo no se está basando fuertemente en atributos que puedan tener un sesgo preocupante de edad o económico, sino que en comportamiento de un cliente. Que las variables se utilicen no es necesariamente algo malo, pues eliminarlas puede hasta aumentar el sesgo., así que no es preocupante que no tengan importancia exactamente 0.

Si bien se determinaron que estas variables potencialmente problemáticas tienen muy poco efecto en el modelo, graficar su contribución puede ayudarnos a entender mejor exactamente como se están usando.

In [40]:
age_index = ebm.term_names_.index("age")
occupation_index = ebm.term_names_.index("occupation")
In [41]:
ebm.explain_global(name="EBM").visualize(age_index)

Primero, hay que recordar que muchas edades en el conjunto de datos no son razonables. Haciendo zoom al rango razonable [0, 120], se notan pequeños efectos de la edad, que no ocurren en otros lados de gráfico. En particular, la edad parece aumentar el riesgo (pero no mucho), exceptuando entre los 45 a 50 años. Este es un sesgo riesgoso, aunque muy pequeño.

In [42]:
ebm.explain_global(name="EBM").visualize(occupation_index)

El trabajo tiene un efecto demasiado minúsculo para distinguirse. Eliminar la variable podría reducir el ruido, dado que no tiene efecto útil.

7. Concluir [1.0 puntos]¶

El mejor modelo solamente tuvo un f1-score de 0.79 en macro promedio. Este resultado es más que satisfactorio, en especial considerando que se empezó con un baseline de solamente 0.29. La resolución del problema fue exitosa, con resultados buenos con un modelo intepretable.

El EDA ayudó a entender la distribución de los datos para crear buenos pipelines para varios modelos, pero los mejores modelos fueron los que realizan codificaciones categóricas por dentro y no requieren de escalamiento, así que la creación de pipelines no fue tan útil. Una gran utilidad que dió el EDA posterior a las transformaciones fue encontrar variables de alta información mutua, que ayudaron a la ingeniería de características para una nueva variable altamente predictiva como el producto de dos variables importantes.

Mejoras que hubiera realizado serían:

  • Considerar interactions=0 para hacer el modelo aún más simple
  • Apagar el bagging de EBM (usando el parámetro outer_bags=1) y utilizar en vez el de imblearn con BalancedBaggingClassifier, para mejorar el desempeño con el desbalance. El problema es que eso puede hacer perder la incertidumbre sobre los puntajes, pero se pueden unir los clasificadores a uno interpretable con interpret.glassbox.merge_ebm.
  • Eliminar características que se determinaron que son irrelevantes en el análisis a posteriori, por ejemplo, usando RFE.
  • Estudiar mejor el experimento de optuna con más hiperparámetros para realizar una mejor guianza en la búsqueda

Utilicé este proyecto como experiencia práctica para aprender de Optuna y más sobre usos prácticos de EBM, que es un modelo que conocía teóricamente, pero no había logrado a usar para intepretaciones. También aproveché de probar la técnica de validación adversarial, que he escuchado pero no implementado con EBM. Me hubiera gustado tener más tiempo para probar más técnicas de optimización de hiperparámetros y de rebalanceo de datos, pero tengo otros trabajos que debo priorizar :(

Created in deepnote.com Created in Deepnote